import os
os.environ["OMP_NUM_THREADS"] = '8'
os.environ["OMP_PROC_BIND"] = "spread" 

import re
import matplotlib.pyplot as plt
import networkx as nx

import numpy as np
np.set_printoptions(precision=7)
import math
import random,time

from projectq.backends import CommandPrinter,Simulator
from projectq.cengines import (MainEngine,DummyEngine,LocalOptimizer,
                               DecompositionRuleSet, BasicEngine, LastEngineException)
import projectq.setups.decompositions
from projectq.setups import restrictedgateset

from projectq.meta import Control
# from projectq.ops import H,Measure,All,Rz,Ry,CNOT,Allocate,CX,Deallocate,CZ
from projectq.ops import *

ANGLE_PRECISION = 0.015

class GatePrinter(BasicEngine):
    """

    """
    def __init__(self, accept_input=False, default_measure=False,
                 in_place=False):

        BasicEngine.__init__(self)
        self.gateList = []

    def is_available(self, cmd):
        """
        Specialized implementation of is_available: Returns True if the
        CommandPrinter is the last engine (since it can print any command).

        Args:
            cmd (Command): Command of which to check availability (all
                Commands can be printed).
        Returns:
            availability (bool): True, unless the next engine cannot handle
                the Command (if there is a next engine).
        """
        try:
            return BasicEngine.is_available(self, cmd)
        except LastEngineException:
            return True

    def _print_cmd(self, cmd):
        """
        Print a command or, if the command is a measurement instruction and
        the CommandPrinter is the last engine in the engine pipeline: Query
        the user for the measurement result (if accept_input = True) / Set
        the result to 0 (if it's False).

        Args:
            cmd (Command): Command to print.
        """
        if(cmd.gate != Allocate and cmd.gate != Deallocate):
            self.gateList.append(str(cmd))
            

    def receive(self, command_list):
        """
        Receive a list of commands from the previous engine, print the
        commands, and then send them on to the next engine.

        Args:
            command_list (list<Command>): List of Commands to print (and
                potentially send on to the next engine).
        """

        for cmd in command_list:
            
            if not cmd.gate == FlushGate():
                self._print_cmd(cmd)
            # (try to) send on
            if not self.is_last_engine:
                self.send([cmd])


class SqrtYGate(BasicGate):
    """ Square-root Y gate class """
    @property
    def matrix(self):
        return 0.5 * np.matrix([[1+1j, -1-1j], [1+1j, 1+1j]])

    def tex_str(self):
        return r'$\sqrt{Y}$'

    def __str__(self):
        return "SqrtY"

SqrtY = SqrtYGate()

class GateOperation(object):
    def __init__(self, gate, target, angle=None):
        self.gate = gate
        self.target = target
        self.angle = angle

    def __eq__(self, value):
        return (
            self.gate == value.gate
            and self.target == value.target
            and self.angle == value.angle
            )

class CircuitDag(object):
    def __init__(self, nQubit, operations):
        self.nQubit = nQubit
        self.dag = self.OpList2Dag(operations)

    def OpList2Dag(self, operations):
        dag = nx.DiGraph()
        borderI = []
        for i in range(self.nQubit):
            borderI.append(GateOperation('borderI', (i,)))
            dag.add_node(-i-1, op=borderI[i])

        for i, op in enumerate(operations):
            dag.add_node(i, op=op)

        lenOps = len(operations)
        for i,op in enumerate(borderI):
            dag.add_node(i+lenOps, op=op)

        endIndex = lenOps+len(borderI)
        for i in range(-len(borderI), endIndex):
            for t in dag.node[i]['op'].target:
                for j in range(i+1, endIndex):
                    if(t in dag.node[j]['op'].target):
                        dag.add_edge(i, j, qid = t)
                        break

        return dag

    def Dag2OpList(self, dag):
        topoOperationList = nx.topological_sort(dag)
        opList = []
        for i in topoOperationList:
            op = dag.node[i]['op']
            if(op.gate != 'borderI'):
                opList.append(op)

        return opList

    def RemoveCircuitNodes(self, dag, nodes):
        for node in nodes:
            inEdges = list(dag.in_edges(node))
            outEdges = list(dag.out_edges(node))
            pattern = r'%d%d'%(len(inEdges),len(outEdges))
            edgePattern = ('11','12','21','22')
            if(len(inEdges) == 0 or len(outEdges) == 0):
                pass
            elif(pattern == edgePattern[3]):
                inNode = (inEdges[0][0], inEdges[1][0])
                outNode = (outEdges[0][1], outEdges[1][1])
                if(dag[inNode[0]][node]['qid'] != dag[node][outNode[0]]['qid']):
                    outNode = outNode[::-1]
                dag.add_edge(inNode[0], outNode[0],
                                    qid = dag[inNode[0]][node]['qid'])
                dag.add_edge(inNode[1], outNode[1],
                                    qid = dag[inNode[1]][node]['qid'])
            elif(pattern == edgePattern[2]):
                inNode = (inEdges[0][0], inEdges[1][0])
                outNode = outEdges[0][1]
                dag.add_edge(inNode[0], outNode,
                                    qid = dag[inNode[0]][node]['qid'])
                dag.add_edge(inNode[1], outNode, 
                                    qid = dag[inNode[1]][node]['qid'])
            elif(pattern == edgePattern[1]):
                inNode = inEdges[0][0]
                outNode = (outEdges[0][1], outEdges[1][1])
                dag.add_edge(inNode, outNode[0], 
                                    qid = dag[node][outNode[0]]['qid'])
                dag.add_edge(inNode, outNode[1],
                                    qid = dag[node][outNode[1]]['qid'])
            elif(pattern == edgePattern[0]):
                inNode = inEdges[0][0]
                dag.add_edge(inNode, outEdges[0][1], 
                                    qid = dag[inNode][node]['qid'])
            
            dag.remove_node(node)

    def TansformGateNodeToEdge(self, dag, node, edge, nodeAttr):
        '''
        Transform unneeded gate node to a single qubit gate node and insert to a edge
        '''
        if(dag.has_node(node) and dag.has_edge(*edge)):
            if(node not in edge):
                self.RemoveCircuitNodes(dag, [node])
                dag.add_node(node, op = nodeAttr)
                dag.add_edge(edge[0], node, qid = nodeAttr.target[0])
                dag.add_edge(node, edge[1], qid = nodeAttr.target[0])
                dag.remove_edge(*edge)

def RestoreMap(mapper, state):
    state = np.array(state)
    stateOut = np.zeros(len(state), dtype = state.dtype)
    mapperBin = []
    for i in mapper:
        mapperBin.append(1<<i)

    mapl = []
    for i in range(len(state)):
        tempSum = 0
        iBin = bin(i)[:1:-1]
        for j,Bin in enumerate(iBin):
            if(Bin == '1'):
                tempSum += mapperBin[j]
        stateOut[i] = state[tempSum]
        mapl.append(tempSum)

    return stateOut,mapl

def InnerProduct(vec1, vec2):
    return np.conj(np.array(vec1).reshape(1,(len(vec1)))).dot(
                    np.array(vec2).reshape((len(vec2),1))
                )[0][0]


def load_circuit(filename):
    """Load the quantum circuit from a file.
    Note:
    You do not need to change this function.
    """
    with open(filename, 'r') as circuit_file:
        circuit = circuit_file.readlines()
        circuit_file.close()
    return circuit

def GetFidelity(CMat, CxList):
    # initlize values
    OneGateDirt = {'I':0.9995,
                'X':0.9992,
                'Y':0.9993,
                'Z':0.9994,
                'H':0.9991,
                'T':0.9997,
                'Rx':0.9993,
                'X/2':0.9993,
                'Ry':0.9995,
                'Y/2':0.9995,
                'Rz':0.9998,
                'Z/2':0.9998}
    FidelityMatrix = []
    
    TimeZero = 10 ** 6
    BitNum = len(CMat)
    OperationDepth = len(CMat[0])
    
    # initlize time S matrix
    TimeSM = [[20] * OperationDepth for i in range(BitNum)]
    
    for i in range(len(CxList)):
        CxBitNum = len(CxList[i]) - 1
        Depth = int(CxList[i][0])
        for j in range(1, CxBitNum + 1):
            BitPos = int(CxList[i][j])
            TimeSM[BitPos][Depth] = CxBitNum * 20

    # Get useful time S and circuit matrix to calculate fidelity
    SumTime = [0] * BitNum
    DepthList = [0] * BitNum
    CountInsert = [0] * BitNum

    for i in range(len(CxList)):
        CxBitNum = len(CxList[i]) - 1
        for j in range(1, CxBitNum + 1):
            BitPos = int(CxList[i][j])
            DepthList[BitPos] = int(CxList[i][0]) + CountInsert[BitPos]
            SumTime[BitPos] = np.sum(TimeSM[BitPos][0:DepthList[BitPos]])
        MaxTime = np.max([SumTime[int(CxList[i][j])] for j in range(1, CxBitNum + 1)])
                
        for j in range(1, CxBitNum + 1):
            BitPos = int(CxList[i][j])
            Depth = DepthList[BitPos]
            for k in range(int((MaxTime - SumTime[BitPos] )/20)):
                TimeSM[BitPos].insert(Depth+k, 20)
                CMat[BitPos].insert(Depth+k, 'I')
                CountInsert[BitPos] += 1

    # calculate fidelity matrix 
    # fidelity = 1
    SumTime = [0] * BitNum
    for i in range(BitNum):
        FidelityMatrix.append([])
        for j in range(len(CMat[i])):
            SumTime[i] += TimeSM[i][j]
            BetaS = math.exp(-SumTime[i]/TimeZero)
            if (CMat[i][j] == 'Tg'):
                FideRS = 1
            elif(CMat[i][j] == 'I'):
                FideRS = 1
            elif(CMat[i][j] == 'Cg' ):
                if (TimeSM[i][j] > 80):
                    raise ValueError('Error, controled gate   is more than four,'
                          'fidelity is less zero')
                else:
                    FideRS = (1-i*0.0001) * BetaS * (1 - 
                                   (10**(TimeSM[i][j]/20)) * 0.000056)

            elif(CMat[i][j] in OneGateDirt):
                FideRS = (1-i*0.0001) * BetaS * OneGateDirt[CMat[i][j]]
            else:
                FideRS = (1-i*0.0001) * BetaS * 0.9995
            # fidelity = fidelity * FideRS
            FidelityMatrix[i].append(FideRS)
   
    # print(CMat)
    # print(TimeSM)
    for i in range(BitNum):
        for j in range(int((max(SumTime) - SumTime[i])/20)):
            TimeSM[i].append(20)
            CMat[i].append('I')
            SumTime[i] += 20
            BetaS = math.exp(-SumTime[i]/TimeZero)
            FideRS = (1-i*0.0001) * BetaS * 0.9995
            # fidelity = fidelity * FideRS
            FidelityMatrix[i].append(FideRS)
    
    # print(CMat)
    # print(TimeSM)
    fidelity = 1
    for FideR in FidelityMatrix:
        for FideRS in FideR:
           fidelity = fidelity * FideRS
      
    return fidelity


def calculate_circuit_fidelity(qubits, circuit):
    """The function used to calculate the fidelity of given circuit.
    Args:
    qubits(int):
    The qubit number of these two quantum circuit.
    circuit(str list):
    The quantum circuit you want to calculate the fidelity with.
    Return:
    fidelity(float):
    The estimated fidelity of the given quantum circuit.
    """
    
    nQubit, operations = GetFormatCircuit(circuit)
    cMat,cxList = GetCircuitMatrix(nQubit, operations)
    
    fidelity = GetFidelity(cMat, cxList)

    return fidelity



def circuit_simplifier(eng, initial_circuit):
    """The function you need to modify only.
     Args:
       eng:
         The MainEngine of HiQ and you do not need to change.
       initial_circuit(str list):
         The initial quantum circuit is loaded from txt files such as
         'circuit1.txt'. An example of given quantum circuit return from
         function load_circuit(filenam) is:
           ['Allocate | Qureg[0]\n', 'Rx(1.570796326795) | Qureg[0]\n',
           'Deallocate | Qureg[0]']
     Return:
       simplified_circuit(str list):
         The simplified quantum circuit using your algorithm in the same
         format as initial_circuit.
     Note:
       The amplitude result simulated from simplified_circuit should be the
       same with that of initial_circuit.
    """

    simplified_circuit = ''
    nQubit, operations = GetFormatCircuit(initial_circuit)
    # # ####################################
    # print('Before Length:')
    # print(len(operations))
    # # wf1 = GetWaveFunction(nQubit, operations)


    operations = ProjectqOptimizer(nQubit, operations)

    circuitDag = CircuitDag(nQubit, operations)

    HardamardMergeRule(circuitDag)
    CxRule(circuitDag)
    CommutationMergeRule(circuitDag)
    CxRule(circuitDag)
    PhasePolynomialsRule(circuitDag, 0)
    CxRule(circuitDag)
    CommutationMergeRule(circuitDag)
    HardamardMergeRule(circuitDag)
    PhasePolynomialsRule(circuitDag, -1)
    HardamardMergeRule(circuitDag)
    CxRule(circuitDag)
    # FinalRule(circuitDag)

    operations = circuitDag.Dag2OpList(circuitDag.dag)
    #CX change to CZ
    temp = []
    for op in operations:
        if(op.gate == 'CX'):
            temp.append(GateOperation('H', (op.target[1],)))
            op.gate = 'CZ'
            temp.append(op)
            temp.append(GateOperation('H', (op.target[1],)))
        else:
            temp.append(op)
    operations = temp
    circuitDag = CircuitDag(nQubit, operations)
    HardamardMergeRule(circuitDag)
    FinalRule(circuitDag)
    operations = circuitDag.Dag2OpList(circuitDag.dag)

    # cMat, cxList = GetCircuitMatrix(nQubit, operations)
    # PlotCircuitMatrix(cMat, cxList)
    # print('Depth after:')
    # print(len(cMat[0]))
    simplified_circuit = OpList2StrCircuit(operations)
    # ###########################################
    # wf2 = GetWaveFunction(nQubit, operations)
    # ii = InnerProduct(wf1,wf2)
    # print('innerProduct:%s'%str(np.abs(ii)))

    return simplified_circuit



def GetRandomCircuit(qubits, NumberGates):
    from random import randint
    from math import pi
    GateList = ['H','CX','Rz','Z','X','Y','T']
    AnglesList = [4*pi, 3*pi, 2*pi, pi, pi/2,pi/4]
    BitsList = []

    f = open('Randomcircuit.txt', 'w')
    for i in range(qubits):
        BitsList.append(i)
        temp = 'Allocate | Qureg['+ str(i) +']\n'
        f.write(temp)

    
    for i in range(NumberGates):
        gate = GateList[randint(0, 2)]
        if(gate == 'Rz'):
            angle = AnglesList[randint(0, 5)]
            qbit = BitsList[randint(0, qubits-1)]
            temp = 'Rz('+str(angle)+') | Qureg[' + str(qbit) + ']\n'
            
        elif(gate == 'CX'):
            while True:
                qbit1 = BitsList[randint(0, qubits-1)]
                qbit2 = BitsList[randint(0, qubits-1)]
                if(qbit1 != qbit2):
                    break

            temp = 'CX | ( Qureg[' + str(qbit1) + '], Qureg[' + str(qbit2) + '] )\n'
        else:
            qbit = BitsList[randint(0, qubits-1)]
            temp = gate + ' | Qureg[' + str(qbit) + ']\n'
        f.write(temp)

    for i in range(qubits):
        temp = 'Deallocate | Qureg['+ str(qubits-1-i) + ']\n'
        f.write(temp)
    f.close()


def PlotCircuitMatrix(circuitMat, cxList, circuitRange = None):
    r = len(circuitMat)
    if(circuitRange==None):
        s = [0, len(circuitMat[0])]
    else:
        s = circuitRange
    dpi = 4000/(s[1]-s[0])
    dpi = 20
    size=(s[1]-s[0], r)
    # size = (100,r)
    plt.figure(figsize = size, dpi = dpi)

    for i in range(r):
        plt.plot([s[0]-1, s[1]],[-i, -i],'k-',alpha=0.8)

    for i in range(r):
        for j in range(s[0], s[1]):
            if(circuitMat[i][j] != 'Cg' and circuitMat[i][j] != 'Tg' and circuitMat[i][j] != 'I'):
                plt.plot([j],[-i],'s',markersize=25,color='#cfcfcf',alpha=1)
                plt.text(j, -i, circuitMat[i][j], fontsize = 14, color='k',
                            ha='center', va='center')

    # sort cxList by x position
    temp = [cx[0] for cx in cxList]
    sortIndex = sorted(enumerate(temp), key=lambda x:x[1])
    temp = []
    for i in sortIndex:
        temp.append(cxList[i[0]])
    cxList = temp

    tempX = 0
    shiftX = 0
    for cx in cxList:
        if(cx[0] != tempX):
            tempX = cx[0]
            shiftX = tempX

        if(cx[0] >= s[0] and  cx[0] < s[1]):
            plt.plot([shiftX, shiftX],[-cx[1], -cx[-1]],'k-',alpha=1)
            for y in cx[1:-1:1]:
                plt.plot([shiftX],[-y],'ro',markersize=5,alpha=1)
            plt.plot([shiftX],[-cx[-1]],'ro',markersize=20,alpha=1)
        shiftX += 0.08

    plt.ylim((-r,1))
    plt.show()


def CircuitFileFilter(filename):
    circuit = []
    temp = load_circuit(filename)
    for s in temp:
        if('#' in s):
            pass
        elif('|' in s):
            circuit.append(s)

    return circuit


def GetFormatCircuit(opStrList):
    '''
    Format the circuit string list to regular operation list\n
    return:\n
        nQubit: int\n
        operations: list of object of [GateOperation] 
    '''

    circuit = []
    temp = opStrList
    for s in temp:
        if('#' in s):
            pass
        elif('|' in s):
            circuit.append(s)

    operations = []
    excludeGateDict = {str(Allocate):1, str(Deallocate):2, str(Measure):3}
    allQid = set()
    try:
        for op in circuit:
            temp = op.split('|')
            gate = re.findall(r'(.*)\((.*)\)', temp[0])
            if(len(gate)): gate = list(gate[0])
            else: gate = [temp[0].strip()]

            for i in range(len(gate)):
                gate[i] = gate[i].strip()

            target = re.findall(r'\[(.*?)\]', temp[1])
            target = tuple([int(i) for i in target])
            
            bitCut = False
            # for i in target:
            #     if i>15:
            #         bitCut = True
            if(not bitCut):
                allQid.update(target)
                if(len(gate)==1):
                    if(not (gate[0] in excludeGateDict)):
                        operations.append(GateOperation(gate[0], target))
                else:
                    if(gate[0] == 'C'):
                        operations.append(GateOperation(temp[0].strip(), target))
                    else:
                        operations.append(GateOperation(gate[0], target, float(gate[1])))
    except Exception:
        raise ValueError(
            "Txt to Gate Operation Transform Error: command: %s; gate: %s"%(op,str(gate)))
    
    nQubit = max(allQid)+1
    return nQubit, operations

def OpList2StrCircuit(operations):
    circuit2 = []
    for op in operations:
        if(op.angle != None):
            circuit2.append(op.gate + r'(%s)'%str(op.angle) + '|' + r'Qureg[%d]'%(op.target[0]))
        elif(len(op.target) == 2):
            circuit2.append(op.gate + '|' + r'(Qureg[%d], Qureg[%d])'%(op.target[0], op.target[1]))
        elif(len(op.target) == 1):
            circuit2.append(op.gate + '|' + r'Qureg[%d]'%(op.target[0]))

    return circuit2

def GetCircuitMatrix(nQubit, operations):
    circuitMat = [[] for i in range(nQubit)]
    CxList = []
        
    for op in operations:
        if(len(op.target) == 1): 
            circuitMat[op.target[0]].append(op.gate)
        else:
            maxLen = 0
            maxQubitId = 0
            for qubitId in op.target:
                if(len(circuitMat[qubitId]) > maxLen):
                    maxLen = len(circuitMat[qubitId])
                    maxQubitId = qubitId

            for qubitId in op.target:
                if(qubitId != maxQubitId):
                    circuitMat[qubitId].extend(['I']*(maxLen-len(circuitMat[qubitId])))
                if(qubitId == op.target[-1]):
                    circuitMat[qubitId].append('Tg')
                else:
                    circuitMat[qubitId].append('Cg')

            CxList.append([maxLen]+list(op.target))

        
    MaxDepth = np.max([len(circuitMat[i]) for i in range(nQubit)])
    for i in range(nQubit):
        circuitMat[i].extend(['I']*(MaxDepth - len(circuitMat[i])))
        
    return  circuitMat, CxList


def PlotGraph(G):
    pos=nx.spring_layout(G) 

    nx.draw_networkx_nodes(G,pos) 
    nx.draw_networkx_labels(G,pos) 
    nx.draw_networkx_edges(G,pos) 
    plt.show()

def PlotSubDagCircuit(circuitDag, nodes):
    dag = circuitDag.dag
    nQubit = circuitDag.nQubit
    subG =  nx.DiGraph(nx.freeze(dag.subgraph(nodes)))
    operations = circuitDag.Dag2OpList(subG)
    cMat, cxList = GetCircuitMatrix(nQubit, operations)
    PlotCircuitMatrix(cMat, cxList)

def GetWaveFunction(nQubit, operations):
    eng = MainEngine(Simulator(), [])

    Qureg = eng.allocate_qureg(nQubit)
    All(H) | Qureg
    try:
        for op in operations:
            target = tuple([Qureg[i] for i in op.target])
            gate = eval(op.gate)
            if(op.angle != None):
                gate(op.angle) | target
            else:
                gate | target
    except Exception:
        raise ValueError("Not imported Gate: %s"%op.gate)
    eng.flush()

    mapper,wf = eng.backend.cheat()
    wf0,_ = RestoreMap(mapper.values(), wf)

    All(Measure) | Qureg

    return np.array(wf0)

class MapGraph(object):

    def __init__(self, mapper):
        self.Graph, self.MaxQuregLength = self.undirected_graph(mapper)

    def undirected_graph(self, mapper):
        G = nx.path_graph(0)
        lines = mapper.split('\n')
        n = int(lines[0])
        lines = lines[1:]
        for line in lines:
            a, b = line.split(',')
            nx.add_path(G, [int(a) - 1, int(b) - 1])
        return G, n

    def find_shortest_path(self, sour, targ):
        list_path = nx.shortest_path(self.Graph, source=sour, target=targ)
        # print('%d节点到%d节点最短路径: ' % (sour, targ), list_path)
        return list_path

    def TransformCnot(self, c, t):
        path = self.find_shortest_path(c, t)
        n = len(path)
        if (n == 2):
            return [[c, t]]
        ctPairList = []
        cnotList = []

        for i in range(n - 1):
            ctPairList.append([path[i], path[i + 1]])
        cnotList += ctPairList
        cnotList += ctPairList[-2::-1]
        cnotList += ctPairList[1:]
        cnotList += ctPairList[-2:0:-1]
        return cnotList

def map_generation(qubit):
    map = str(qubit)
    for i in range(qubit):
        map = map + '\n' + str(i) + ',' + str(i+1)
    return map


def LastProjectqOptimizer(nQubit, operations):
    engList = restrictedgateset.get_engine_list((Rx,Ry,Rz,H),(CX,))
    eng = MainEngine(GatePrinter(), engList)

    mapper = map_generation(nQubit+1)
    mapG = MapGraph(mapper)

    def CnotTranslation(control, target):
        paths = mapG.TransformCnot(control, target)
        for path in paths:
            c = path[0]-1
            t = path[1]-1
            CX | (Qureg[c], Qureg[t])

    Qureg = eng.allocate_qureg(nQubit)
    try:
        for op in operations:
            gate = eval(op.gate)
            target = tuple([Qureg[i] for i in op.target])
            if (op.angle != None):
                if(op.gate == 'Rx'):
                    H | target
                    Rz(op.angle) | target
                    H | target
                elif(op.gate == 'Ry'):
                    Rz(np.pi*7/2) |target
                    H | target
                    Rz(op.angle) |target
                    H |target
                    Rz(np.pi/2) | target
                elif(op.gate == 'Rz'):
                    angle = op.angle
                    if(angle > 1): angle = 4*np.pi - angle
                    if(angle > ANGLE_PRECISION):
                        gate(op.angle) | target
                else:
                    gate(op.angle) | target
            else:
                if(op.gate == 'X'):
                    H | target
                    Rz(np.pi) | target
                    H | target
                elif(op.gate == 'T'):
                    Rz(np.pi/4) | target
                elif (op.gate == 'SqrtX'):
                    H | target
                    Rz(np.pi / 2) | target
                    H | target
                elif (op.gate == 'SqrtY'):
                    Rz(np.pi) | target
                    H | target
                elif(op.gate == 'CZ' or op.gate == 'C(Z)'):
                    H | target[1]
                    (control, target) = op.target
                    CnotTranslation(control+1, target+1)
                    H | target[1]
                elif(op.gate == 'CX'):
                    (control, target) = op.target
                    CnotTranslation(control+1, target+1)
                else:
                    gate | target
    except Exception:
        raise ValueError("Not imported Gate: %s"%op.gate)
    eng.flush()

    _, newOps = GetFormatCircuit(eng.backend.gateList)
    All(Measure) | Qureg

    return newOps


def ProjectqOptimizer(nQubit, operations):
    engList = restrictedgateset.get_engine_list((Rx,Ry,Rz,H),(CX,))
    eng = MainEngine(GatePrinter(), engList)

    Qureg = eng.allocate_qureg(nQubit)
    try:
        for op in operations:
            gate = eval(op.gate)
            target = tuple([Qureg[i] for i in op.target])
            if (op.angle != None):
                if(op.gate == 'Rx'):
                    H | target
                    Rz(op.angle) | target
                    H | target
                elif(op.gate == 'Ry'):
                    Rz(np.pi*7/2) |target
                    H | target
                    Rz(op.angle) |target
                    H |target
                    Rz(np.pi/2) | target
                elif(op.gate == 'Rz'):
                    angle = op.angle
                    if(angle > 1): angle = 4*np.pi - angle
                    if(angle > ANGLE_PRECISION):
                        gate(op.angle) | target
                else:
                    gate(op.angle) | target
            else:
                if(op.gate == 'X'):
                    H | target
                    Rz(np.pi) | target
                    H | target
                elif(op.gate == 'T'):
                    Rz(np.pi/4) | target
                elif (op.gate == 'SqrtX'):
                    H | target
                    Rz(np.pi / 2) | target
                    H | target
                elif (op.gate == 'SqrtY'):
                    Rz(np.pi) | target
                    H | target
                elif(op.gate == 'CZ' or op.gate == 'C(Z)'):
                    H | target[1]
                    CX | target
                    H | target[1]
                else:
                    gate | target
    except Exception:
        raise ValueError("Not imported Gate: %s"%op.gate)
    eng.flush()

    _, newOps = GetFormatCircuit(eng.backend.gateList)
    All(Measure) | Qureg

    return newOps

def HardamardMergeRule(circuitDag):
    dag = circuitDag.dag
    allNodes = list(dag.nodes)
    for i in allNodes:
        if(dag.has_node(i)):
            if(dag.node[i]['op'].gate == 'CX'):
                ## H1-H2-CX-H1-H2
                neighbors = list(dag.predecessors(i)) + list(dag.successors(i))
                if(len(neighbors) == 4):
                    nH = 0
                    for n in neighbors:
                        if(dag.node[n]['op'].gate == 'H'): nH += 1
                        else: break
                    if(nH == 4):
                        circuitDag.RemoveCircuitNodes(dag, neighbors)
                        dag.node[i]['op'].target = dag.node[i]['op'].target[::-1]
                
                ## H2-S2-CX-Sdag2-H2
                target = dag.node[i]['op'].target[-1]
                psNode = []
                for node in list(dag.predecessors(i)):
                    if(dag[node][i]['qid'] == target): psNode.append(node)
                for node in list(dag.successors(i)):
                    if(dag[i][node]['qid'] == target): psNode.append(node)
                if(len(psNode)==2):
                    if(dag.node[psNode[0]]['op'].gate == 'Rz'
                        and dag.node[psNode[1]]['op'].gate == 'Rz'):
                        angle1 = dag.node[psNode[0]]['op'].angle 
                        angle1 = angle1 if angle1 < 2*np.pi else angle1 - 4*np.pi
                        angle2 = dag.node[psNode[1]]['op'].angle
                        angle2 = angle2 if angle2 < 2*np.pi else angle2 - 4*np.pi
                        if(abs(angle1+angle2) < ANGLE_PRECISION
                            and abs(abs(angle1)-np.pi/2) < ANGLE_PRECISION ):
                            nodeH = list(dag.predecessors(psNode[0])) + list(dag.successors(psNode[1]))
                            if(len(nodeH) == 2):
                                if(dag.node[nodeH[0]]['op'].gate == 'H'
                                    and dag.node[nodeH[1]]['op'].gate == 'H'):
                                    circuitDag.RemoveCircuitNodes(dag, nodeH)
                                    dag.node[psNode[0]]['op'].angle = \
                                        4*np.pi - angle1 if angle1>0 else -angle1
                                    dag.node[psNode[1]]['op'].angle = \
                                        4*np.pi - angle2 if angle2>0 else -angle2

            #H-H
            elif(dag.node[i]['op'].gate == 'H'):
                removenodes = list()
                go_on = True
                currentNode = i
                while(go_on):
                    node0 =  list(dag.successors(currentNode))[0]
                    if(dag.node[node0]['op'].gate == 'H'):
                        removenodes.append(node0)
                        currentNode = node0
                    else:
                        go_on = False
                if(len(removenodes) % 2 == 1):
                    removenodes.append(i)
                circuitDag.RemoveCircuitNodes(dag, removenodes)


def CommutationMergeRule(circuitDag):
    dag = circuitDag.dag

    allNodes = list(dag.nodes)
    # Rz-CX-Rz-CX-Rz
    for node in allNodes:
        if dag.has_node(node):
            if dag.node[node]['op'].gate == 'Rz':
                target_qubit = dag.node[node]['op'].target[0]
                neighbors_pre, neighbors_suc = list(dag.predecessors(node)), list(dag.successors(node))
                if len(neighbors_pre) and len(neighbors_suc):
                    node_pre, node_suc = neighbors_pre[0], neighbors_suc[0]
                    if dag.node[node_pre]['op'].gate == 'CX'\
                        and dag.node[node_pre]['op'] == dag.node[node_suc]['op']\
                        and dag.node[node_pre]['op'].target[1] == target_qubit:
                        neighbors_pre, neighbors_suc = list(dag.predecessors(node_pre)), \
                                                       list(dag.successors(node_suc))
                        if len(neighbors_pre) == 2 and len(neighbors_suc) == 2:
                            if(dag[neighbors_pre[0]][node_pre]['qid'] == target_qubit):
                                neighbors_pre = neighbors_pre[::-1]
                            if(dag[node_suc][neighbors_suc[0]]['qid'] == target_qubit):
                                neighbors_suc = neighbors_suc[::-1]
                            node_pre_pre = neighbors_pre[1]
                            node_suc_suc = neighbors_suc[1]
                            if dag.node[node_pre_pre]['op'].gate == 'Rz'\
                                    and dag.node[node_suc_suc]['op'].gate == 'Rz':
                                node_CX_suc = list(dag.successors(node_pre))
                                if len(node_CX_suc) == 2:
                                    if node_CX_suc[0] == node:
                                        next_node = node_CX_suc[1]
                                    else:
                                        next_node = node_CX_suc[0]
                                    go_on = True
                                    while go_on:
                                        if next_node == node_suc:
                                            dag.node[node_pre_pre]['op'].angle\
                                                = (dag.node[node_pre_pre]['op'].angle
                                                   + dag.node[node_suc_suc]['op'].angle) % (4 * np.pi)
                                            circuitDag.RemoveCircuitNodes(dag, [node_suc_suc])
                                            go_on = False
                                        elif dag.node[next_node]['op'].gate == 'CX':
                                            go_on = False
                                        elif dag.node[next_node]['op'].gate == 'Rz':
                                            next_node = list(dag.successors(next_node))[0]
                                            go_on = True

    allNodes = list(dag.nodes)
    # Rz-H-Cx-H-Rz
    for node in allNodes:
        if dag.has_node(node):
            if dag.node[node]['op'].gate == 'CX':
                target_qubit = dag.node[node]['op'].target[1]
                neighbors_pre, neighbors_suc = list(dag.predecessors(node)), list(dag.successors(node))
                if len(neighbors_pre) ==2 and len(neighbors_suc) == 2:
                    if(dag[neighbors_pre[0]][node]['qid'] == target_qubit):
                        neighbors_pre = neighbors_pre[::-1]
                    if(dag[node][neighbors_suc[0]]['qid'] == target_qubit):
                        neighbors_suc = neighbors_suc[::-1]
                    node_pre = neighbors_pre[1]
                    node_suc = neighbors_suc[1]
                    if dag.node[node_pre]['op'].gate == 'H' \
                            and dag.node[node_suc]['op'].gate == 'H':
                        neighbors_pre, neighbors_suc = list(dag.predecessors(node_pre)), \
                                                       list(dag.successors(node_suc))
                        if len(neighbors_pre) == 1 and len(neighbors_suc) == 1:
                            node_pre_pre, node_suc_suc = neighbors_pre[0], neighbors_suc[0]
                            if dag.node[node_pre_pre]['op'].gate == 'Rz' \
                                    and dag.node[node_suc_suc]['op'].gate == 'Rz':
                                angle_pre, angle_suc = dag.node[node_pre_pre]['op'].angle, \
                                                       dag.node[node_suc_suc]['op'].angle
                                circuitDag.RemoveCircuitNodes(dag, [node_suc_suc])
                                dag.node[node_pre_pre]['op'].angle = (angle_pre + angle_suc) % (4*np.pi)
    
    allNodes = list(dag.nodes)
    # Cx-H-Cx-H-Cx
    for node in allNodes:
        if dag.has_node(node):
            if dag.node[node]['op'].gate == 'CX':
                neighbors = list(dag.predecessors(node)) + list(dag.successors(node))
                if len(neighbors) == 4:
                    node_pre, node_suc = neighbors[1], neighbors[2]
                    if dag.node[node_pre]['op'].gate == 'H' \
                            and dag.node[node_suc]['op'].gate == 'H':
                        neighbors_pre, neighbors_suc = list(dag.predecessors(node_pre)), \
                                                       list(dag.successors(node_suc))
                        if len(neighbors_pre) == 1 and len(neighbors_suc) == 1:
                            node_pre_pre, node_suc_suc = neighbors_pre[0], neighbors_suc[0]
                            if dag.node[node_pre_pre]['op'] == dag.node[node_suc_suc]['op'] \
                                    and dag.node[node_pre_pre]['op'].gate == 'CX' \
                                    and len(dag.node[node_pre_pre]['op'].target) == 2:
                                if dag.node[node_pre_pre]['op'].target[1] == dag.node[node_pre]['op'].target[0] \
                                        and dag.node[node_pre_pre]['op'].target[0] != dag.node[node]['op'].target[1]:
                                    edges_CX_suc = list(dag.out_edges(node_pre_pre))
                                    control_CX = dag.node[node_pre_pre]['op'].target[0]
                                    if len(edges_CX_suc) == 2:
                                        removeNodes = []
                                        for edge in edges_CX_suc:
                                            if dag[edge[0]][edge[1]]['qid'] == control_CX:
                                                next_node = edge[1]
                                                go_on = True
                                                while go_on:
                                                    if next_node == node_suc_suc:
                                                        removeNodes = [node_pre_pre, node_suc_suc]
                                                        go_on = False
                                                        break
                                                    elif dag.node[next_node]['op'].gate == 'CX':
                                                        go_on = False
                                                    elif dag.node[next_node]['op'].gate == 'Rz':
                                                        next_node = list(dag.successors(next_node))[0]
                                                        go_on = True   
                                        circuitDag.RemoveCircuitNodes(dag, removeNodes)


def CxRule(circuitDag):
    dag = circuitDag.dag
    allNodes = list(dag.nodes)
    for i in allNodes:
        if(dag.has_node(i)):

            # Rz-Rz
            if(dag.node[i]['op'].gate == 'Rz'):
                removenodes = list()
                go_on = True
                currentNode = i
                startNode = i
                angles = dag.node[i]['op'].angle
                while(go_on):
                    nodes0 =  list(dag.successors(currentNode))
                    for node0 in nodes0:
                        if(dag.node[node0]['op'].target[0] == dag.node[startNode]['op'].target[0]):
                            if(dag.node[node0]['op'].gate == 'Rz'):
                                removenodes.append(node0)
                                angles += dag.node[node0]['op'].angle
                                currentNode = node0
                                go_on = True
                                break
                            elif(dag.node[node0]['op'].gate == 'CX'):
                                currentNode = node0
                                break
                            else:
                                go_on = False
                                break
                        else:
                            go_on = False
                            

                dag.node[i]['op'].angle = angles % (4*np.pi) 
                circuitDag.RemoveCircuitNodes(dag, removenodes)
                # CaXb-Rza-CaXb
                prenode = list(dag.predecessors(i))[0]
                sucnode =  list(dag.successors(i))[0]
                if(dag.node[prenode]['op'].gate == 'CX'
                    and dag.node[prenode]['op'] == dag.node[sucnode]['op']
                    and dag.node[prenode]['op'].target[0] == dag.node[i]['op'].target[0]
                    and dag.has_successor(prenode, sucnode)):
                    circuitDag.RemoveCircuitNodes(dag, [prenode, sucnode])
            
            # 
            elif(dag.node[i]['op'].gate == 'CX'):
                removenodes = list()
                target1 = dag.node[i]['op'].target
                prenodes = list(dag.predecessors(i)) 
                sucnodes =  list(dag.successors(i))
                                
                # CaXb-CaXb-
                if (len(prenodes) == 1
                        and dag.node[prenodes[0]]['op'] == dag.node[i]['op']):
                    removenodes = [prenodes[0], i]
                # -CaXb-CaXb
                elif (len(sucnodes) == 1
                        and dag.node[sucnodes[0]]['op'] == dag.node[i]['op']):
                    removenodes = [i, sucnodes[0]]

                # 
                elif(len(prenodes) == 2 and len(sucnodes) == 2):
                    for node1 in prenodes:
                        for node2 in sucnodes:
                            if(dag.node[node1]['op'].gate == 'CX'
                            and dag.node[node1]['op'] == dag.node[node2]['op']):
                                # CaXb-CaXn-CaXb
                                if(dag.node[node1]['op'].target[0] == target1[0] 
                                and dag.has_successor(node1, node2)):
                                    removenodes = [node1, node2]
                                # CaXb-CnXb-CaXb
                                elif(dag.node[node1]['op'].target[1] == target1[1] 
                                and (dag.has_successor(node1, node2) 
                                or dag.node[node1]['op'].gate == 'Rz')):
                                    removenodes = [node1, node2]
                circuitDag.RemoveCircuitNodes(dag, removenodes)
  

def FinalRule(circuitDag):
    dag = circuitDag.dag
    allNodes = list(dag.nodes)
    for i in allNodes:
        if(dag.has_node(i)):
            # CX-X1-CX = X1-X2 
            if(dag.node[i]['op'].gate == 'X'):
                node1 = list(dag.predecessors(i))[0]
                node2 = list(dag.successors(i))[0]
                if(dag.node[node1]['op'].gate == 'CX' and 
                        dag.node[node1]['op'] == dag.node[node2]['op'] and 
                        dag.node[node1]['op'].target[0] == dag.node[i]['op'].target[0] and
                        dag.has_successor(node1,node2)):
                    op = GateOperation('X', (dag.node[node1]['op'].target[1],))
                    inNodes = list(dag.predecessors(node1))
                    circuitDag.RemoveCircuitNodes(dag, [node1])
                    for n in inNodes:
                        if(dag.has_successor(n, node2)):
                            dag.add_node(node1, op = op)
                            dag.add_edge(n, node1, qid = op.target[0])
                            dag.add_edge(node1, node2, qid = op.target[0])
                            dag.remove_edge(n, node2)
                            break
                    circuitDag.RemoveCircuitNodes(dag, [node2])

            if(dag.node[i]['op'].gate == 'Rz'):
                angle1 = dag.node[i]['op'].angle
                angle1 = angle1 if angle1 < 2*np.pi else angle1 - 4*np.pi
                if(abs(abs(angle1)-np.pi) < ANGLE_PRECISION):
                    node1 = list(dag.predecessors(i))[0]
                    node2 =  list(dag.successors(i))[0]
                    # Rz(pi) replace by Z
                    dag.node[i]['op'].gate = 'Z'
                    dag.node[i]['op'].angle = None
                    # H-Rz(pi)-H = X
                    if (dag.node[node1]['op'].gate == 'H' and
                            dag.node[node2]['op'].gate == 'H'):
                        dag.node[i]['op'].gate = 'X'
                        dag.node[i]['op'].angle = None
                        circuitDag.RemoveCircuitNodes(dag, [node1, node2])
                    # CX-Rz(pi)2-CX    
                    elif(dag.node[node1]['op'].gate == 'CX' and 
                            dag.node[node1]['op'] == dag.node[node2]['op'] and 
                            dag.node[node1]['op'].target[1] == dag.node[i]['op'].target[0] and
                            dag.has_successor(node1, node2)):
                        op = GateOperation('Z', (dag.node[node1]['op'].target[0],))
                        inNodes = list(dag.predecessors(node1))
                        circuitDag.RemoveCircuitNodes(dag, [node1])
                        for n in inNodes:
                            if(dag.has_successor(n, node2)):
                                dag.add_node(node1, op = op)
                                dag.add_edge(n, node1, qid = op.target[0])
                                dag.add_edge(node1, node2, qid = op.target[0])
                                dag.remove_edge(n, node2)
                                break
                        circuitDag.RemoveCircuitNodes(dag, [node2])

                # elif((abs(abs(angle1)-np.pi/2) < ANGLE_PRECISION)):
                #     dag.node[i]['op'].gate = 'S'
                #     dag.node[i]['op'].angle = None

                elif((abs(abs(angle1)-np.pi/4) < ANGLE_PRECISION)):
                    dag.node[i]['op'].gate = 'T'
                    dag.node[i]['op'].angle = None

def GetPhasePolynomialSubGraph(dag, origin_node):
    def CxSearch(dag, Cx_node, current_qubit, direction):
        operate_qubit = dag.node[Cx_node]['op'].target
        if(direction == 'L'):
            in_or_out_edges = list(dag.in_edges(Cx_node)) 
        elif(direction == 'R'): 
            in_or_out_edges = list(dag.out_edges(Cx_node)) 

        if len(in_or_out_edges) == 2:
            for edge in in_or_out_edges:
                if dag[edge[0]][edge[1]]['qid'] == current_qubit:
                    if(direction == 'L'): next_node = edge[0]
                    elif(direction == 'R'): next_node = edge[1]
                    return next_node, set(operate_qubit)
        elif len(in_or_out_edges) == 1:
            edge = in_or_out_edges[0]
            if(direction == 'L'): next_node = edge[0]
            elif(direction == 'R'): next_node = edge[1]
            return next_node, set(operate_qubit)
        return None
    borderList = ('H', 'borderI')
    sub_circuit = {}
    xxxxx = []
    sub_Cx = []
    
    qubit_list = set()
    origin_qubit = dag.node[origin_node]['op'].target
    current_qubit = origin_qubit[1]
    qubit_list.add(current_qubit)
    N = 0

    while len(qubit_list) > N:
        N = len(qubit_list)
        row_circuit = []
        if dag.node[origin_node]['op'].gate == 'CX':
            row_circuit.append(origin_node)
            sub_Cx.append(origin_node)
            new_qubit = []
            next_Cx = []
            origin_tc = set(dag.node[origin_node]['op'].target)

            left_node, _ = CxSearch(dag, origin_node, current_qubit, 'L')
            left_gate = dag.node[left_node]['op'].gate
            while left_gate not in borderList:
                row_circuit.append(left_node)
                if left_gate == 'CX':
                    sub_Cx.append(left_node)
                    temp = left_node
                    left_node, operate_qubit = CxSearch(dag, left_node, current_qubit, 'L')
                    left_gate = dag.node[left_node]['op'].gate
                    if len(operate_qubit-origin_tc) != 0:
                        next_qubit = operate_qubit - qubit_list
                        for q in next_qubit:
                            new_qubit.append(q)
                            next_Cx.append(temp)
                else:
                    node_list = list(dag.predecessors(left_node))
                    left_node = node_list[0]
                    left_gate = dag.node[left_node]['op'].gate
            row_circuit.append(left_node)

            right_node, _ = CxSearch(dag, origin_node, current_qubit,'R')
            right_gate = dag.node[right_node]['op'].gate
            while right_gate not in borderList:
                row_circuit.append(right_node)
                if right_gate == 'CX':
                    sub_Cx.append(right_node)
                    temp = right_node
                    right_node, operate_qubit = CxSearch(dag, right_node, current_qubit, 'R')
                    right_gate = dag.node[right_node]['op'].gate
                    if len(operate_qubit - origin_tc) != 0:
                        next_qubit = operate_qubit - qubit_list
                        for q in next_qubit:
                            new_qubit.append(q)
                            next_Cx.append(temp)
                else:
                    node_list = list(dag.successors(right_node))
                    right_node = node_list[0]
                    right_gate = dag.node[right_node]['op'].gate
            row_circuit.append(right_node)

            sub_circuit[current_qubit] = row_circuit
            if len(new_qubit) == 0:
                current_qubit = origin_qubit[0]
                if current_qubit not in qubit_list:
                    qubit_list.add(current_qubit)
                    origin_node = sub_Cx[0]
                else:
                    N = len(qubit_list)
            else:
                current_qubit = new_qubit[0]
                origin_node = next_Cx[0]
                qubit_list.add(current_qubit)

    delet_Cx = {'C':[], 'T':[]}
    for node in sub_Cx:
        count = 0
        isinCT = [False, False]
        for qubit, row_circuit in sub_circuit.items():
            if node in row_circuit:
                count = count + 1
                if(dag.node[node]['op'].target[0] == qubit):
                    isinCT[0] = True
                else: isinCT[1] = True
        if count != 2:
            if(isinCT[0]):  delet_Cx['C'].append(node)
            if(isinCT[1]):  delet_Cx['T'].append(node)
 
    return sub_circuit, delet_Cx, sub_Cx


def PhasePolynomialGateMerge(circuitDag, subCircuitNodeList, delet_Cx):
    dag =  circuitDag.dag
    outCx_Cin = delet_Cx['C']
    outCx_Tin = delet_Cx['T']
    outCx = outCx_Cin + outCx_Tin
    border = ('H', 'borderI')
    leftBorder = []
    nodes = []
    maxQubit = 0
    for l in subCircuitNodeList.values():
        nodes.extend(l)
        for n in l:
            if(dag.node[n]['op'].gate in border):
                leftBorder.append(n)
                break

    subG =  nx.DiGraph(nx.freeze(dag.subgraph(nodes)))
    topoNodes = list(nx.topological_sort(subG))
    phaseDict = {}
    def addPhase(phase, nodeORedge):
        if(phase not in phaseDict):
            phaseDict[phase] = {'Rz':[], 'edge':[]}
        if(type(nodeORedge) != tuple):
            phaseDict[phase]['Rz'].append(nodeORedge)
        else:phaseDict[phase]['edge'].append(nodeORedge)
            
    for n in leftBorder: # initial: left border out edge add PhasePolynomials
        outNode = list(subG.successors(n))
        qid = subG[n][outNode[0]]['qid']
        phase = 2**qid
        if(maxQubit < qid): maxQubit = qid
        subG[n][outNode[0]]['P'] = phase
        if(subG.node[outNode[0]]['op'].gate != 'Rz'):
            addPhase(phase, (n,outNode[0]))

    maxQubit += 1
    # handle outCx whose target bit is in subG
    for cx in outCx_Tin:
        outNode = list(subG.successors(cx))
        for node in outNode:
            if(subG.node[node]['op'].gate not in border):
                subG[cx][node]['P'] = 2**maxQubit

    for n in topoNodes:
        if(subG.node[n]['op'].gate not in border):
            inNode = list(subG.predecessors(n))
            outNode = list(subG.successors(n))
            control = subG.node[n]['op'].target[0]
            if(subG[inNode[0]][n]['qid'] != control):
                inNode = inNode[::-1]
            if(subG[n][outNode[0]]['qid'] != control):
                outNode = outNode[::-1]
            # handle outCx whose control bit is in subG
            if(n in outCx_Cin):
                for node in outNode:
                    if('P' in subG[inNode[0]][n]):
                        subG[n][node]['P'] = subG[inNode[0]][n]['P']
                    else:
                        pass
            elif(subG.node[n]['op'].gate == 'CX'
                    and n not in outCx_Tin):
                pattern = (len(inNode), len(outNode))
                if(pattern[1] == 2):
                    if(pattern[0] == 2):
                        inPhaseCT = (subG[inNode[0]][n]['P'], subG[inNode[1]][n]['P'])
                    elif(pattern[0] == 1): 
                        inPhaseCT = subG[inNode[0]][n]['P']
                        if(subG.node[n]['op'] != subG.node[inNode[0]]['op']):
                            inPhaseCT = inPhaseCT[::-1] # two sequential inverse CXs
                    subG[n][outNode[0]]['P'] = inPhaseCT[0]
                    if( np.log2(inPhaseCT[0]+1e-5) >= maxQubit and np.log2(inPhaseCT[1]+1e-5) >= maxQubit):
                        subG[n][outNode[1]]['P'] = inPhaseCT[0] + inPhaseCT[1]
                    else:
                        subG[n][outNode[1]]['P'] = inPhaseCT[0] ^ inPhaseCT[1]
                    for node in outNode:
                        if(subG.node[node]['op'].gate != 'Rz'):
                            phase = subG[n][node]['P']
                            if(phase < 2**maxQubit): addPhase(phase, (n, node))
                elif(pattern[1] == 1):
                    if(pattern[0] == 2):
                        inPhaseCT = (subG[inNode[0]][n]['P'], subG[inNode[1]][n]['P'])
                    elif(pattern[0] == 1):
                        inPhaseCT = subG[inNode[0]][n]['P']
                        if(subG.node[n]['op'] != subG.node[inNode[0]]['op']):
                            inPhaseCT = inPhaseCT[::-1] # two sequential inverse CXs
                    if( np.log2(inPhaseCT[0]+1e-5) >= maxQubit and np.log2(inPhaseCT[1]+1e-5) >= maxQubit):
                        subG[n][outNode[0]]['P'] = (inPhaseCT[0], inPhaseCT[0] + inPhaseCT[1])
                    else:
                        subG[n][outNode[0]]['P'] = (inPhaseCT[0], inPhaseCT[0] ^ inPhaseCT[1])
                    for phase in subG[n][outNode[0]]['P']:
                        if(phase < 2**maxQubit): addPhase(phase, (n, outNode[0]))
            elif(subG.node[n]['op'].gate == 'Rz'):
                phase = subG[inNode[0]][n]['P']
                subG[n][outNode[0]]['P'] = phase
                if(phase < 2**maxQubit): addPhase(phase, n)

    #Check and merge CXs 
    for n in topoNodes:
        if(subG.node[n]['op'].gate == 'Rz'):
            target = subG.node[n]['op'].target[0]
            inNode = list(subG.predecessors(n))
            outNode = list(subG.successors(n))
            if(len(inNode) and len(outNode)
                    and subG.node[inNode[0]]['op'].gate == 'CX'
                    and subG.node[outNode[0]]['op'] == subG.node[inNode[0]]['op']
                    and subG.has_successor(inNode[0], outNode[0])
                    and inNode[0] not in outCx
                    and outNode[0] not in outCx
                    and subG.node[inNode[0]]['op'].target[1] == target
                    and subG.node[outNode[0]]['op'].target[1] == target):
                for key, value in phaseDict.items():
                    if(n in value['Rz']):
                        if(len(value['Rz']) > 1):
                            phaseDict[key]['Rz'].remove(n)
                            angle = (dag.node[value['Rz'][0]]['op'].angle 
                                        + dag.node[n]['op'].angle)
                            angle = angle % (4*np.pi)
                            dag.node[value['Rz'][0]]['op'].angle = angle
                            circuitDag.RemoveCircuitNodes(dag, [n, inNode[0], outNode[0]])
                        elif(len(value['edge']) > 0):
                            for edge in value['edge']:
                                if(dag.node[edge[0]]['op'].gate != 'CX'
                                    or dag.node[edge[1]]['op'].gate != 'CX'):
                                    nodeAttr = dag.node[n]['op']
                                    nodeAttr.target = (dag[edge[0]][edge[1]]['qid'],)
                                    circuitDag.TansformGateNodeToEdge(dag, n, edge, nodeAttr)
                                    circuitDag.RemoveCircuitNodes(dag, [inNode[0], outNode[0]])
                                    phaseDict[key]['edge'].remove(edge)
                                    break

    # PlotSubDagCircuit(dag, set(nodes) - set(outCx))
    #Merge Rz
    RzMergeList = []
    for value in phaseDict.values():
        if(len(value['Rz']) > 1):
            RzMergeList.append(value['Rz'])
    for Rzs in RzMergeList:
        angle = 0
        for Rz in Rzs:
            angle += dag.node[Rz]['op'].angle
        angle = angle % (4*np.pi)
        dag.node[Rzs[0]]['op'].angle = angle
        circuitDag.RemoveCircuitNodes(dag, Rzs[1:])


def PhasePolynomialsRule(circuitDag, direction):
    dag = circuitDag.dag
    allNodes = list(dag.nodes)
    CXList = set()
    for node in allNodes:
        if dag.node[node]['op'].gate == 'CX':
            CXList.add(node)

    while len(CXList) > 0:
        # randInt = random.randint(0,len(CXList)-1)
        origin_node = list(CXList)[direction]

        sub_circuit, delet_Cx, sub_Cx = GetPhasePolynomialSubGraph(dag, origin_node)
        PhasePolynomialGateMerge(circuitDag, sub_circuit, delet_Cx)
        CXList = CXList - set(sub_Cx)
 
if __name__ == "__main__":
    filename = 'circuit4.txt'

    t1 = time.time()
    circuit = CircuitFileFilter(filename)
    fidelityBefore = calculate_circuit_fidelity(0, circuit)
    print('F_before: %s'%str(fidelityBefore))

    simpcircuit = circuit_simplifier(None, circuit)
    print('Simplifier Length:')
    print(len(simpcircuit))

    fidelityAfter = calculate_circuit_fidelity(0, simpcircuit)
    print('F_After: %s'%str(fidelityAfter))
    

    print(time.time()-t1)